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Abstract Adding sand grains at a single site in Abelian sandpile models produces beautiful but 
complex patterns. We study the effect of sink sites on such patterns. Sinks change the scaling of the 
diameter of the pattern with the number N of sand grains added. For example, in two dimensions, 
in presence of a sink site, the diameter of the pattern grows as -^Z {N/ log A) for large A, whereas it 
grows as VA if there are no sink sites. In presence of a line of sink sites, this rate reduces to A^/^. 
We determine the growth rates for these sink geometries along with the case when there are two lines 
of sink sites forming a wedge, and its generalization to higher dimensions. We characterize one such 
asymptotic patterns on the two-dimensional F- lattice with a single source adjacent to a line of sink 
sites, in terms of position of different spatial features in the pattern. For this lattice, we also provide 
an exact characterization of the pattern with two sources, when the line joining them is along one of 
the axes. 
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1 Introduction 

It is well known that beautiful and complex patterns can be generated by deterministic evolution 
of systems under simple local rules, e.g. in the game of life HI, and Turing patterns [5]. Growing 
sandpiles on a flat table with boundaries by adding particles at a constant rate gives rise to singular 
structures like ridges in the stationary state, which have attracted much attention recently [3lll]- In the 
Abelian sandpile model, growing sandpiles produce richer and hence more interesting patterns. This 
model is inspired by real sandpile dynamics, but has different rules of evolution. The steady state of 
sandpile models with slow driving, and presence of a boundary has been studied much in the context 
of self-organized criticality [5]- The Abelian sandpile model, with particles added at one site, on an 
infinite lattice have the very interesting property of proportionate growth [6]. This is a well-known 
feature of biological growth in animals, where different parts of the growing animal grow at roughly 
the same rate. Our interest in studying growing sandpiles comes from it being the prototypical model of 
proportionate growth. Most of the other growth models studied in physics literature, such as diffusion- 
limited aggregation, or surface deposition do not show this property, and the growth is confined to 
some active outer region, and the inner structures once formed are frozen in, and do not evolve further 
in time [7]. 

In [6], we studied growing sandpiles in the Abelian model on the F-lattice and the Manhattan 
lattice. These are directed variants of the square lattice, obtained by assigning directions to the bonds, 
as shown in Fig.l. We fomid that for a particular choice of the initial background configuration, the 
patterns formed can be characterized exactly. The special initial configuration is the one in which 
each alternate site of the lattice is occupied, forming a chequerboard pattern. If we add particles at 
the origin, and relax the configuration using the sandpile toppling rules, we generate a fairly complex 
pattern made up of triangles and dart-shaped patches (Fig[2]), that shows proportionate growth. The 
full characterization of this pattern reveals an interesting underlying mathematical structure, which 
seems to deserve further exploration. This is what we do in this paper, by adding sink sites, or multiple 
sources. 

Presence of sink sites changes the pattern in interesting ways. In particular, it changes how different 
spatial lengths in the pattern scale with the number of added grains N . For example, in absence of 
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Fig. 1 Directed square lattices studied in this paper;(a) F-lattice and (b) Manhattan lattice. 



sink sites, the diameter of the pattern grows as for large A'', whereas presence of a single sink 
site next to the site of addition, this changes to a N/ log growth. If there is a line of sink sites 
next to the site of addition the growth rate is N^^^. We also studied the case where the source site 
is at the corner of a wedge-shaped region of wedge angle uj = 7r/2, 37r/2, or 27r, and where the wedge 
boundaries are absorbing. (The last case corresponds to the source next to an infinite half line.) For 
the single point source, determination of different distances in the pattern required a solution of the 
Laplace equation on a discrete Riemann surface of two sheets. Interestingly, for these wedge angles, 
we still have to solve the discrete Laplace equation, but the structure of the Riemann surface changes, 
e.g. from two-sheets to three-sheets for lu ^ it, and five-sheets for lu = 2it. We characterize the patterns 
in terms of the solution of the discrete Laplace equation. We also show that the pattern grows as TV", 
with a = 2lo/{tt + Auj). 

We also study the effect of having multiple sites of addition on the pattern. For multiple sources, 
the pattern of small patches near each source is not substantially different from a single-source pattern, 
but some rearrangement occurs in the larger outer patches. Two patches may sometimes join into one, 
or conversely, a patch may break up into two. But the number of patches undergoing such changes 
is finite. However, the sizes of all patches are affected by the presence of other sources, and we show 
how these changes can be calculated exactly for the asymptotic pattern. Spatial patterns in sandpile 
models were first discussed by Liu et al [8]. The asymptotic shape of the boundaries of sandpile 
patterns produced by adding grains at single site on different periodic backgrounds was discussed 
in [9]. Borgne [TOj obtained bounds on rate of growth of these boundaries and later these bounds 
are improved by Fey et al |llj and Levine et al [12j . The first detailed analysis of different periodic 



Fig. 2 Stable configuration for tlie ASM, obtained by adding 5 x 10'' grains at one site, on tiie F-lattice of 
FigUJa) witli initial chequerboard configuration. Color code: red=0, yellow=l. Apparent orange regions in the 
picture represent patches with chequerboard configuration. (Details can be seen in the online version using 
zoom in.) 

structures found in the patterns are carried out by Ostojic in |13| . Other special configurations in the 
ASM models, like the identity [T0 l [T4 l [T5] . or the stable state produced from special unstable states 
also show complex internal self-similar structures [8|, which also share common features with the 
patterns studied here. There are other models, which are related to the Abelian sandpilc model, e.g. 
the internal Diffusion-Limited Aggregation (DLA), Eulerian walkers (also called rotor-router model), 
and the infinitely-divisible sandpilc, which also show similar structure. For the Internal DLA, Gravncr 
and Quastel showed that the asymptotic shape of the growth pattern is related to the classical Stefan 
problem in hydrodynamics, and determined the exact radius of the pattern with a single point source 
[16j . Levine and Peres have studied patterns with multiple sources in these models recently, and proved 
the existence of a limit shapeflT]. 

This paper is organized as follows. After defining the model in Section 2, we discuss scaling of the 
diameter of the patterns with N. We first consider in Section 3 the pattern in the presence of a line of 
sink sites. In Section 4, this analysis is extended to other sink geometries: two intersecting line sinks 
in two dimensions and two or three intersecting plane sinks in three dimensions. The case of a single 
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sink site is a bit different from others, and is discussed separately in Section 5. The remaining sections 
are devoted to a detailed characterization of some of these patterns. In Section 6 we give a summary 
of our earlier work on characterization of single source pattern, and use it to characterize the pattern 
in presence of a line sink. In Section 7, we discuss the case when there arc two sources present. Section 
8 contains a summary and some concluding remarks. 

2 Definition of the model 

We consider the Abelian sandpilc model on the F-lattice (Fig[T]2). This is a square lattice with directed 
bonds such that each site has two inward and two outward arrows. A different assignment of arrow 
directions, that gives us the Manhattan lattice is shown in FigllJ. The asymptotic pattern formed by 
growing sandpilc on the Manhattan lattice is the same as on the F-lattice [6]. We shall discuss here 
only the F-latticc, but the discussion is equally applicable to the Manhattan lattice. 

Define a position vector on the lattice, R = {x, y). In the Abelian sandpilc model, a height variable 
z (R), called the number of grains on the site, is assigned to each site R. In a stable configuration all 
sites have height z (R) < 2. The system is driven by adding grains at a single site and if this addition 
makes the system unstable it relaxes by the toppling rule: each unstable site transfers one grain each 
in the direction of its outward arrows. We start with an initial configuration in which z (R) = 1, for 
sites with {x + y) = even, and otherwise. For numerical purpose we used a lattice large enough so 
that none of the avalanches reaches the boundary. The result of adding iV = 5 x 10^ grains at the 
origin is shown in FiglS] 

3 Growth of the pattern with line sink 

Consider the pattern formed by adding sand grains at a single site in presence of a line of sink sites. 
Any grain reaching a sink site gets absorbed, and is removed from the system. For simplicity let us 
consider the source site at Ro = (2^07 0) and the sink sites along the y-axis. A picture of the pattern 
produced by adding 14336000 grains at (1, 0) is shown in FigO When N grains have been added, let 
2A (N) be the diameter of the pattern, measured as the height of the smallest rectangle that enclosed 
all sites that have toppled at least once. We want to study how A{N) increases as a function of N. 



Fig. 3 Pattern produced by adding grains at a single site adjacent to a line of sink sites. Color code: red=0 and 
yellow=l. Apparent orange regions in the picture represent patches with chequerboard configuration. (Zoom 
in for details in the online version.) 

As mentioned before, the pattern exhibits proportionate growth. While there is as yet no rigorous 
proof of this important property, we assume this in the following. Then, it is natural to describe the 
pattern in reduced coordinates defined hy ^ ^ x/A and rj ~ y/A. A position vector in this reduced 
coordinate is defined by r = Tl/ A = (^,77). Then in A —>■ 00, the pattern can be characterized by a 
function Ap{r) which gives local excess density of sand grains in the pattern in a small rectangle of 
size 6£,5ri about the point r, with 1/A-^ 5^, Sr] <C 1. 

Let T/i (R) be the number of toppling at site R when the diameter reaches value 2A for the first 
time. Define 

where R' = ([vi^J, [Ari\) with [.tJ being the floor function which gives the largest integer < x. 

From the conservation of sand-grains in the toppling process, it is easy to see that (j) satisfies the 
Poisson equation [6] 

V^(bir) = Ap{r)-^S{r~ro) (2) 

for all r in the right-half plane with ^ > 0, where is the position of the source in reduced coordinates. 
Also as there are no toppling at sink sites, (j) niust satisfy the boundary condition 

(/) (r) = for all r = (0, t]) . (3) 
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A complete information of (r) determines the density function Ap (r) and in turn characterizes the 
asymptotic pattern. 

We can think of cj) as the potential due to a point charge N/A^ at Tq and an arcal charge density 
—Ap (r), in presence of a grounded conducting line along the ?7-axis. This problem can be solved using 
the well-known method of images in electrostatics. Let r' be the image point of r with respect the 
77-axis. Define Ap (r) in the left half plane as 

Ap{r')^-Apir). (4) 

Then the Poisson equation for this new charge configuration is 

V20(r) = Apir) ~^Hr~ r„) + - r'o) . (5) 

As the function Ap (r) is odd under reflection. automatically vanishes along the 77-axis. 
We define Nr as the number of sand grains that remain unabsorbed. Then 

7V. = ^^ziz(a:,y), (6) 

x>0 y 

where Az (x, y) is the change in height variables before and after the system relaxes. Clearly, for large 
A, we can write 

N,.-A^[dTAp{r), (7) 

where dr = d^drj is the infinitesimal area around r = (^,77) and the integration performed over the 
right half-plane H with ^ > 0. We shall use the sign ~ to denote equality up to leading order in A. Since 
Ap (r) is a non-negative bounded function, exactly zero outside a finite region, this integral exists. Let 
its value be C2 and then we have 

Nr ^C2A^. (8) 

Let Na denote the number of grains that are absorbed by the sink sites. Then considering that 
grains can reach sink sites only by toppling at its neighbors we have 

Na^\Y.TA{l,y). (9) 

y 

The factor 1/2 comes from the fact that in F-lattice, only half of the sites on the column x = \ would 
have arrows going out to the sink sites. Then using our scaling ansatz in equation for A large, 
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Hence 

iV, ^A' r dr,^ . (11) 

Now from equation ([5]) the potential (f) can be written as sum of two terms: 4>dipoie due to two point 
charges N/A^ and —N/A^ at Tq = (^o, 0) and its image point rj, = (— ^o, 0) respectively, and the term 
4>rest due to the areal charge density. 

(j) (r) (j)dipole (r) + (j)rest (r) , (12) 

where 

y^(l)dipoie (r) = (r - To) + (r - r;,) , 

V^(j)rest{r)^Ap{r). (13) 

We first consider the case where Ro is finite and = Ro/ A vanishes in the large A limit. Then 4>dipoie 
reduces to a dipole potential, and it diverges near origin. However, (prest (r) is a non-singular function 
for all r. From the solution of dipole potential, it is easy to show that 

(l)dipoie{r,e)K a'^^^, (14) 
r 

for 1 3> r ^ 1/yl, where we have used polar coordinates (r, 0) with 9 being measured with respect 
to the ^-axis. Here A is a numerical constant, which is determined by the property of the asymptotic 
pattern. Then 
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4' (15) 



and the integral in equation (jlip diverges as A/rjmin, where r]min is the cutoff introduced by the lattice. 
Using rjmin = O (1/^) it is easy to show that 

Na ^ CiA\ (16) 

where Ci is a constant. Then using equation ^ and (flBl) and that Na and A^^ add up to N , we get 

CiA^ + C2A^ = N. (17) 

Here we use the symbol = to denote "nearly equal to" . Considering the dominant term in the expression 
for large A, it follows that A increases as N^^^. 
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The above scaling behavior is verified with our numerical data. Let A* (N) be the real positive root 
of the equation p?)) for a given value of integer N. As A takes only integer values on the lattice, an 
estimate of it would be Mint [A*{N)] the integer nearest to A* (N). 

Interestingly, we found that for a choice of C\ = 0.1853 and C2 = 0.528, this estimate gives values 
which differ from the measured value A (N) at most by 1 for all N in the range of 100 to 3 x 10^. 
Clearly more precise estimates of Ci and C2 would be required if we want this to work for larger N. 
Here wc find Eqs.® and (|16p on dimensional counting grounds, and the final Ea. (|17p is then only a 
statement of conservation of sand grains. It is quite remarkable that this scaling analysis gives almost 
the exact value of A{N). The equation has an important feature. It includes "correction to scaling" 
term whereas the usual scaling analysis ignores the sub-leading powers. 

For patterns in the other limit where the source is placed at a distance O (A) such that is non 
zero for A 00, (pdipoie is non-singular along the sink line. Then, clearly Na ^ A^ and as a result 
A{N) - 

4 Generalization to more complex patterns 

The above analysis can be easily generalized to a case with sink sites along two straight lines intersecting 
at an angle uj and a point source inside the wedge. For square lattice, uj — Q, 7r/2, tt, 37r/2 and 2tt are 
most easily constructed, and avoid problems of lines with irrational slopes, or slopes of rational numbers 
with large denominators. The wedge with wedge-angle cj = 7r/2 is obtained by placing sink sites along 
the X and j/-axis and the source site at Ro = (1, 1) in the first quadrant. The pattern with line sink, 
discussed in previous section, correspond to a; = tt. 

For any general lo, corresponding electrostatic problem reduces to determining the potential function 
(/) inside a wedge formed by two intersecting grounded conducting lines. Again the potential has two 
contributions: the potential (ppoint (r) due to a point charge at the source site and the potential (prest (r) 
due to the areal charge density. We first consider the case where the source site is placed at a finite 
distance from the wedge corner such that the distance in reduced coordinate vanishes for A large limit. 
In this limit (firest is non-singular function of r while (f>point diverges close to the origin. A simple 
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calculation of the electrostatic problem gives 

^ , ^, .sina^ 

(l>poMir,9) A—^, (18) 

where a = tt/lu and we have used polar coordinates (r, 6) with the polar angle 6 measured from one of 
the absorbing lines. Again ^ is a constant independent of iV or yl and is a property of the asymptotic 
pattern. Then arguing as before, we get 

Na-CiA^+°' and N, ~ CsA^. (19) 

So the equation analogous to equation p7|) is 

Ciyl^+" + Cavl^ 4= N. (20) 

For wedge angle a; = tt, a = 1, and the above equation reduces to Eq. (fT7|) . For u; ~ 2tt, the value of a 
is 1/2. Again. Eq. (PO)) is in very good agreement with our numerical data. Let A*{N) be the solution 
of Eq. dini) for a given N. Choosing Ci = 0.863408 and C2 = 0.043311, we find that the function 
Nint [A* (N)] differ from the measured values of A at most by 1 for all N in the range from 100 to 
2 X 10^ 

For the problem where the source site is at a distance O (A) from the wedge corner both the 
functions (jirest and (ppoint are nonsingular close to the origin. Then it is easy to show that A (N) grows 
as iVi/2. 

The argument is easily extended to other lattices with different initial height distributions, or to 
higher dimensions. Consider, for example, an Abelian sandpile model defined on the cubic lattice. 
Allowed heights are to 5, and a site topples if the height exceeds 5, and sends one particle to each 
neighbor. The sites are labelled by the Cartesian coordinates {x,y,z), where x,y and z are integers. 
We consider the infinite octant defined by a; > 0, y > 0, z > 0. We start with all heights 4, and add 
sand grains at the site (1,1,1). We assume that the sites on planes x = 0, y ^ and z = are all 
sink sites, and any grain reaching there is lost. We add N grains and determine the diameter of the 
resulting stable pattern. 

We again reduce the potential function in two parts: (fipoint due to a point charge at {1/ A, 1/ A, 1/ A) 
and (jirest due to bulk charge density in presence of three conducting grounded planes. Then simple 
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electrostatic calculation gives that the potential (ppoint is the octapolar potential of the form 

(21) 

where the spherical coordinate is used to denote position. This then implies that the equation deter- 
mining the dependence of yl on is 

CiA^ + C-iA^ 4= N (22) 

Like the other cases, this relation is also confirmed against numerical data. The function Nint [A* {N)] 
with Ci = 0.0159 and C2 = 88 gives almost exact values of A [N). We have checked for N between 
5 X 10^ to 5 X 10*, the difference is at most 1. 

5 A single sink site 

Let the site of addition is the origin and the sink site is placed at Rq. We will show that when To lies 
in a dense patch (color yellow in Figd]), the asymptotic patterns are identical to the one produced in 
absence of a sink. 

The patterns produced for Tq close to 1 with sink sites placed deep inside a dense patch is simple 
to analyze, even for finite but large A. One such pattern is presented in FiglH 

We see that the effect of sink site on the pattern is to produce a depletion pattern centered at the 
sink site. The depletion pattern is a smaller copy of the single source pattern. We define the function 
Azsink (R; N) as the difference between the heights at R in the final stable configuration produced by 
adding N grains at the origin, with and without sink. 

Azs.,nk (R; N) = Az 

source-\-sink 

(R;iV)-Z\z,<,„,,,(R;iV). (23) 

From the figure it is seen that, in this case, Azsink (R; N) is negative of the pattern produced by a 
smaller source, centered at Rq. The number of particles required to produce this smaller pattern is 
exactly the number of particles Na absorbed at the sink site. 

AZs^nk (R; N) = -AZsource (R - Roi A^a) • (24) 

This is immediately seen from the fact that the toppling function Ta (R) satisfies 
TA{Ii')-2TA{^) = Az 

source-\-sink 

(R; N) - NSr^o + NaS-R^R^ , (25) 

R'Gr(R) 
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Fig. 4 Pattern produced by adding 224000 grains at the origin with a sink site at (400, 0) inside a patch 
of density 1 (color yellow ). Color code red= and yellow= 1. The apparent orange regions correspond to 
chequerboard height distribution. (Zoom in for details in the online version.) 




Fig. 5 Pattern produced by adding 224000 grains at origin with a sink site placed at (360, 140) inside a low 
density patch. Color code red= and yellow= 1. The apparent orange regions correspond to the chequerboard 
height distribution. (Details can be seen in the online version using zoom in.) 
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where /^(R) is the set of two neighbors which transfer grains to the site R under topphng. Let 
Tsource(R] N) be the number of topphng at R, when we add N particles at the origin. Since Eq. 
(|25|) is a hnear equation, it fohows that a solution of this equation is 

Ta (R) = Tsource (R; N) — Tgource (R — Ro; ^a) , (26) 

This is a valid solution of our problem, if the corresponding heights in the final configuration with sink 
are all non-negative. This happens when the region with nonzero Azsink is confined within the dense 
patch. 

The number Na can be determined from the requirement that number of toppling at the sink site 
is zero. The potential function for a single source problem diverges as (47r) ^ logr near the source [6]. 
Considering the ultra violate cutoff due to lattice, Tsource (R, N) at R = can be approximated by 
(47r)^^ A^logiV in leading orders in N . Then at R = Ro, Tsource (R — Ro; No) is approximately equal 
to {Any^Na logiVa whereas Tsource (Ro; N) « N(l)source{'ro), where 4> source (r) is the potential function 
for the problem without a sink. Then from equation p6|) we have 

^Na log Na ~ NcPsource {^o) ■ (27) 

Ait 

For large N , this implies that Na — AiKpsource (tq) N/ log A^. In numerical measurement it is found that 
for a change of N from 224000 to 896000, Na\ogN/N changes by less than 7% which is consistent 
with the above scaling relation. For large N, for a sink at a fixed reduced coordinate Tq, the relative 
size of the defect produced by the sink decreases as (log N)^^^^ . Hence asymptotically, the fractional 
area of the defect region will decrease to zero, if sink position Tq is in a dense patch. 

When the sink site is inside a light patch, the subtraction procedure of equation gives positive 
heights, and no longer gives the correct solution. However it is observed for patches in the outer layer 
where patches are large, effect of sink site is confined within neighboring dense patches (Figl5|) and 
rest of the pattern in the asymptotic limit remains unaffected. 

The pattern where the source and sink sites are adjacent to each other appears to be very similar 
to the one produced in absence of sink site. This is easy to see. The Poisson equation analogous to 
equation ([5|) for this problem is 

V20(r) = Ap{r) - ^Mr) + ^S{r - r„), (28) 
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where Na is the number of grains absorbed in the sink site at r,,. In an electrostatic analogy, as discussed 
earlier, (j) can be considered as the potential due to a distributed charge of density —Ap (tq) and two 
point charges of strength N/A^ and —Na/A'^ placed at origin and at respectively. It is easy to see 
that the dominant contribution in the potential is the monopolc term with net charge {N — Na)/ A^. 
The contribution due to other terms decreases as 1/A for large A, and the asymptotic pattern is the 
same as without a sink, with N — Na particles added. 

The number of particles absorbed Na is determined by the condition that the number of toppling 
at (1,0) ( the sink position) is zero. The potential produced by the arcal charge density at (1,0) and 
(0,0) is nearly the same. The number of toppling at (1,0) if we add Na particles at the sink site 
is approximately (Air)'^ Na log Na- Now, from the solution of the discrete Laplacian, the number of 
toppling produced at (1, 0) due to N particles added at (0, 0) is approximately (47r)~^ (A^ log — CN) 
with C being an undetermined constant. Equating these two, we get 

Na\ogNa = N log N -CN (29) 

The above relation is verified with numerical data in figure [H] We find that (A^ logiV — Na log Na)/N 
asymptotically approaches a value C = 2.155 with the difference from the asymptotic value decreasing 
as N~^/^. As the asymptotic pattern is the same as produced by adding {N — Na) grains at the origin 
without a sink, we have N - NaC^ A^. Then, using the numerical value for C we get 

{N - A^) log{N -A^) = N log N - 2.1557V. (30) 

Simplification of this equation for large N, shows that A grows as y/ N / log N with N . 

For finite A, the leading correction to (p (r) comes from the dipolc term in the potential. Presence 
of this term breaks the reflection symmetry of the pattern about the origin. The relative contribution 
of the dipolc potential compared to the monopolc term decays as log A/ A. A measure of the bilateral 
asymmetry is the difference of boundary distances on two opposite sides of the source. This difference 
is plotted in Figl?! where i?i and i?2 arc boundary distances measured along the positive and negative 
X axis with a sink site placed at (1,0). The difference (i?2 — -Ri) is found to fit to 1.22 log (i?2 + 0.5). 
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Fig. 6 Dependence of number of absorbed grains Na on number of grains added A*" at the origin with a sink 
site at (1, 0). 
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Fig. 7 Bilateral asymmetry due to the presence of a sink site in Figl4l 
6 Characterization of the pattern with Une sink 



We start by recalling the characterization of single source pattern [6]. As discussed in Section 3, the 
asymptotic patterns can be characterized by the function Ap (r) in rcscaled coordinate. The single 
source pattern on F-lattice with chequerboard background is made of union of distinct regions, called 
"patches", where inside each patch Ap{r) is constant and takes only two possible values. 1/2 in a 
dense patch and in a light patch [5]. 
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The potential function (r) for the single source problem follows the Poisson equation 

V20(r) = Ap{v) - ^6{v). (31) 

The condition that determines 4> (r) is the requirement that inside each patch of constant density, it is 
a quadratic function of ^ and -q [5]. Let us write 

(j) (r) = + 2h^'q + brf + + eij + f (32) 

where a, h, b, d, e and / are constants inside a patch and a + b = Ap/2 corresponding to the patch. 
Then each patch is characterized by these parameters. Continuity of (r) and its derivatives along the 
boundary between two adjacent patches imposes linear relations among the corresponding parameters. 
These linear equations can be solved on the connectivity graph of patches which forms a square lattice 
on two sheeted Riemann surface [B]. 

The pattern with line sink (FiglSj retained two important properties present in the single source 
pattern. These are: the asymptotic pattern is made of union of two types of patches of excess density 
1/2 and and the separating boundaries of patches are straight lines of slope 0, ±1 or oo. However the 
adjacency graph is changed significantly and this changes the sizes of patches as well. In this section 
we show how to explicitly determine the potential function on this adjacency graph. 

The adjacency graph of the patches is given in FiglHl This representation of the graph is easier to 
see by taking l/?*'^ transformation of the pattern and then joining the neighboring patches by straight 
lines (Fig[H]). Each vertex in the graph is connected to four neighbors except the vertices corresponding 
to the patches next to the absorbing line. These have coordination number 3. Also the vertex at the 
center corresponding to the exterior of the pattern is connected to seven neighbors. 

Let us write the quadratic potential function in a patch P having excess density 1 /2 as 

'^p W = ^('Tip + 1)C^ + ^'^p^^+ ^(1 - 'mp)if +dpi + 6^?? + /p, (33) 

where the parameters m, n, d, e and / take constant values within a patch. Similarly for the lighter 
patches P' 

(0 = Irn^, ie - Tf) + \n^,^i^ + d^,^ + e^,r^ + f^, . (34) 
Using the continuity of (j) (r) and its derivatives along the common boundaries between neighboring 
patches it has been shown that for single source pattern without sink sites m and n take integer values 
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Fig. 8 transformation of the pattern in Fig|3] Two adjoining patches are connected by drawing a straight 

line. 




0,-4 0,-3 0,-2 0,-1 0,0 0,1 0,2 0,3 0,4 



Fig. 9 Adjacency graph of the patches in pattern in FiglS] 

[6]. Same arguments also applies for this problem and (m,n) arc the coordinates of the patches in the 
adjacency graph in Figl^l These coordinates are shown next to some of the vertices. 

There are two different patches corresponding to same set of (771,71) values. Infact, as in the single 
source pattern the adjacency graph forms a square lattice on two sheeted Riemann surface, the same is 
formed for this pattern but on a three sheeted Riemann surface. The pattern covers half of the surface 
with (711, 71) being the Cartesian coordinates on the surface. 

Define function D (m, n) = d (m, n) + ie (?7i, 71) on this lattice. Using the matching conditions along 
the common boundaries between neighboring patches it can be shown that d and e satisfy discrete 
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Cauchy-Riemann condition [6] 



d (m + 1, n + 1) — d (m, n) — e (m, n + 1) — e (m + 1, n) , 



e (m + 1, 71 + 1) — e (m, n) = d {m + l,n) — d (m, rt + 1) , 



(35) 



and then the function D follows discrete Laplaces equation 



D{m + i,n + j) - 4£'(m, n) = 0, 

i=±l j=±l 



(36) 



on this adjacency graph. Let us define M = m + in and z = ^ + t]. Then, as argued before, close to 
the origin the potential diverges as 1/r (equation (jl4p ). Then the corresponding complex potential 
<P (z) ^ 1/z. As M ~ d'^<P/dz^, and D ^ d(l>/dz, it follows that for large |m| + \n\, 



The condition that on the absorbing line (f> (r) must vanish implies that for the vertices with even n 
along the red line in Figl5]e(0, n) vanishes. These vertices correspond to the patches with absorbing line 
as horizontal boundary in FigO Equation (|36p with above constraint and the boundary condition in 
equation p7|) has a unique solution. The normalization of (p is fixed by the requirement that d{l, 0) = —1 
which fixes the diameter of the pattern to be 2 in reduced units. 

Equation [36] is the standard two-dimensional lattice Laplace equation, whose solution is well-known 
when (m, n) G T? jl9| . In our case when the lattice sites form a surface of two sheets, we have not been 
able to find a closed-form formula for D{m,n). However the solution can be determined numerically 
to very good precision by solving it on a finite grid — L < m,n < L with the above conditions imposed 
exactly at the boundary. The calculation is performed with D = 71/^/"^ at the boundary and then the 
solution is normalized to have d (1, 0) = —1. We determined d and e numerically for L ~ 100, 200, 300, 
400 and 500 and extrapolated our results for L ^ oo. 

Comparison of results from this numerical calculation and that obtained by measurements on the 
pattern is presented in Table 1 . We considered four different lengths i?i , i?2 , ^3 and i?4 on the pattern 
and they are shown in Fig llOl Among them, according to the definition of the diameter of the pattern, 
Ri ~ 2 A. We present values of i?2, ^3 and R4 normalized by Ri for different TV. Theoretical values of 
these lengths are determined from asymptotic values of d and e. Comparision of these results shows 
very good agreement among the theoretical and measured values. 



D ~ m2/3. 



(37) 
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Fig. 10 Spatial lengths _Ri, R2, Rz and _R4 tabulated in Table 1. 



N 


896k 


14336A: 


57344ifc 


229376fc 


Theoretical 


R2 
Ri 


0.769 


0.768 


0.770 


0.770 


0.7698 


R3. 

Ri 


0.675 


0.675 


0.667 


0.668 


0.6666 


Ri 
Ri 


0.609 


0.609 


0.617 


0.616 


0.6172 



Table 1 Comparison of different lengths measured directly from the pattern in Fig llOl for increasing values of 
A'', with their theoretical values. 



7 Patterns with two sources 

In this section wc discuss patterns produced by adding N grains each at two sites placed at a distance 
2Aro from each other along the x-axis at Avg and — ylro with r^ = (Co>0). Again the diameter 2A is 
defined as the height of the smallest rectangle enclosing all sites that have toppled atleast once. Two 
limits To close to zero and vq large are trivial: For Tq — > 0, the asymptotic pattern is same as that 
produced by adding grains at a single site. On the otherhand if To > 1, each source produces its own 
pattern, which do not overlap, and the final pattern is simple superposition of the two patterns. 

As noted before, the connectivity graph for single source pattern has square lattice structure on 
a Riemann surface of two sheets [Hj. Then the graph for two non- intersecting single source pattern is 
square lattice on two disjoint Riemann surfaces each consists of two sheets (Fig |13p . Only the vertex 
at the origin represent the exterior of the pattern, which is same for both the single source patterns. It 
has sixteen neighbors and is placed midway between the two Riemann surfaces. For later convenience 
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Fig. 11 Pattern produced by adding iV = 640000 grains each at (-760,0) and (760,0) on F-lattice with 
initial chequerboard distribution of grains and relaxing. This corresponds to ro = 0.95. Color code red=0 and 
yellow=l. (Details can be seen in the online version using zoom in ) 




Fig. 12 Pattern constructed by combining two single source patterns and drawing connecting lines between 
few patches following the connectivity in the pattern in Fig llll 

let us associate the lower Riemann surface to the pattern around — and denote it by F]^. Similarly 
the upper Riemann surface as Fn corresponding to the pattern around Vq. 

For < To < 1, using the Abelian property, we can first topple as if the other source was absent. 
The resulting pattern still has some unstable sites in the region where the patterns overlap. Further 
relaxing these sites transfers these excess grains outward, and changes the dimensions and positions 
of patches: some patches become bigger, some may merge, and some times a patch may beak into two 
disjoint patches. 
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Fig. 13 Representation of adjacency graph of patches for two non-overlapping single source patterns as square 
grid on two Riemann surfaces each of two sheets. Vertices with the same (m, n) coordinates on different sheets 
are represented by different colors. 

The pattern produced with two sources with tq = 0.95 is shown in Fig[TlJ We see that there are 
still only two types of periodic patches, corresponding to /\p(r) values and 1/2, and the slopes of the 
boundaries between patches takes values 0, ±1 or oo. 

The relaxations due to overlaps change the adjacency graph from the case with no overlap. This 
modified adjacency graph is shown in Fig |14l However, for tq just below 1, these changes are few, and 
are listed below. 

i) We note that patches labelled A and A' in FiglT^lhave the same potential function cj). Then, for 
ro just below 1, these patterns can join with each other by a thin strip. This only requires a small 
movement in the boundaries of nearby patches, ( i.e. only a small change in the d and e values of 
nearby patches). Thus, in the adjacency graph, the vertices corresponding to A and A' are collapsed 
into a single vertex A in FigfHl 

ii) Similarly, the vertices corresponding to patches B and B' in Fig[T2]are collapsed into a single 
vertex B in FigfHl 

iii) This divides the region outside the pattern in three parts, O, O' and O". They arc also shown 
in Fig[T3]as separate vertices. 
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^M-L-i ^»Jo,-31 ^^<l,-2) ^M2,-2) ~^0.-l) 



Fig. 14 Connectivity graph for two intersecting single source patterns around two sites of addition placed at 
a distance 2ro from each other. The graph has the structure of square grids on four Riemann sheets except for 
a finite number of vertices indicated by alphabates A, B, O, O' , O" and D shown placed in middle layer. This 
graph remains unchanged for Vo in the range 0.65 to 1.00. 

iv) The patches marked C and C also have the same quadratic form, and the vertical boundary 
between them disappears. However, the patches D and D' are also joined by a thin strip. This horizontal 
strip divides the joined C and C" into two again fFig fT2|) . 

The adjacency of other patches remains unchanged. The adjacency graph of the pattern is shown in 
Fig ll4l Interestingly, this new adjacency graph remains the same for all 0.70 < tq < 1, even though for 
To < 0.85, the sizes of different patches are substantially different. [Compare the pattern for ro = 0.70 
in FigfTSl with pattern for rg = 0.95 in FiglTT]. Shape of the patches near the center of FigfTSl is 
different from that in FigfTTl 

In Figll41 we have have placed the vertices which are formed by merging or dividing patches in 
the overlap region midway between the Riemann sheets corresponding to the two sources. As vq is 
decreased below 0.70, more collisions between growing patches will occur, and the number of vertices 
in this middle region will increase. For any nonzero ro, the number of vertices in the middle layer is 
finite. In the Tq ^ limit vertices from both the surfaces Fl and come together and form a single 
Riemann surface corresponding to a single source pattern around r = 0. For ro small, but greater than 
zero, the outer patches are arranged as in the single-source case, but closer to the sources, one has a 
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Fig. 15 Pattern produced by adding 640000 grains at site (—600,0) and (600,0). Although the pattern is 
significantly different from the one in Fig llll their adjacency graph is same. 

crowded pattern near each source. In the adjacency graph, this corresponds to vertices near the patch 
(0, 0) roughly arranged as on a Ricmann surface of two sheets, while the ones farther from the patch 
(0, 0) remain undisturbed on the 4-sheeted Riemann surface. 

We now characterize the pattern with two sources, and rg > 0.70 in detail by explicitly determining 
the potential function on this adjacency graph. 

The Poisson equation analogous to Eq. ((3T|l for this problem is 

V2</>(r) = Apir) - ^S{r - r„) - ^J(r + r^). (38) 

Let us use the same quadratic form of the potential function given in equation (j33p and equation (|34p . 

Again using the same argument given in [5] it can be shown that m and n are the coordinates of the 
patches in both the adjacency graphs in Fig |13l and Fig 1 141 These coordinates are shown next to each 
vertex. Also, in the region away from the origin, on each sheet, the fimction D{Tn, n) ~ d{m, ?i)+ie(m, n) 
satisfies the discrete Laplace equation 

X! ^ D{m + i,n + j) -4:D{m,n) ^0. (39) 

i=±l j=±l 

Let us define Zq ^ Co + i^o where (^o, rjo) and (— fo, ^Vo) are the coordinates corresponding to Vq and 
—To- Considering that close to Yq and —Yq the potential (r) diverges logarithmically it can be shown 
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Fig. 16 Spatial lengths Ri, R2, R3, R4, and R5 tabulated in Table. 2 
(as done for single source pattern in [6]) that for large \m\ + \n\, 

D{m, n) ~ ± \/l\7 + Zq — , on Fl 

V 27r 4 

= ±^VM-.-„^, onTH (40) 

where yl is a constant independent of N or A. 

Again we determine the solution of equation (|39p numerically to very good precision by solving 
it on finite grid — L < rn, n < L with the conditions Eq. (I40p imposed exactly at the boundary. The 
value of A is determined from a self consistency condition that the diameter of the pattern in reduced 
coordinate is 2 which imposes 2e(— 1, 0) = — 1 corresponding to the vertex A in FigdU We determined 
d and e numerically for L = 100, 200, 300, 400 and 500 and extrapolated our results for L 00. 

Comparison of results from this numerical calculation and that obtained by measurements on the 
pattern in presented in Table 2. We considered five different lengths in the pattern corresponding to 
To = 0.800. These different lengths are drawn in Figlin] and their values rescaled by VN, for the 
patterns with increasing N, is given in Table 2. Theoretical results are obtained using the asymptotic 
values of d and e for large L. The rescaled lengths extrapolated to infinite N limit matches very well 
to the theoretical results. 
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N 


2.5fc 


lOfc 


40fc 


160 A; 


640fc 


Theoretical 


Ri 
viv 


1.84 


1.84 


1.84 


1.83 


1.83 


1.82 


«2 


1.06 


1.07 


1.07 


1.06 


1.05 


1.06 


Rs 

Vn 


0.22 


0.21 


0.20 


0.19 


0.18 


0.18 


Ri 

Vn 


0.18 


0.19 


0.19 


0.18 


0.18 


0.18 


Rs 
\/jV 


0.20 


0.22 


0.21 


0.21 


0.21 


0.21 



Table 2 Comparison of different lengths measured directly from the two source pattern for ro = 0.800 with 
their theoretical values. 

8 Summary 

We have shovifn that exact characterization of the patterns in F-lattice on chequerboard background 
reduces to solving a discrete Laplace equation on the adjacency graph of the pattern. For the single 
source pattern this graph is a square grid on two-sheeted Riemann surface and in presence of a line sink 
it is on a 3-sheeted Riemann surface. This Riemann surface structure occurs for other sink geometries 
also and the number of sheets can be determined from the way (j) diverges near origin. 

If the potential (f>{r) diverges as r~" near the origin, then the corresponding complex function 
<P{z) - z-". Then ■£j<P ~ z^^-". in all the cases we studied above, the patch to which point z belongs 
is characterized by integers (m, n), where -j^^ ^ m + in. Also ^ d + ie. Writing D = d + ie, and 

l + a 

M = m + in, we see that D ^ M 2+a . This then gives the number of Riemann sheets. For example, 
for a wedge angle lo = 2tt, we have a = 1/2. Then D ^ M'^/^, and the Riemann surface would have 5 
sheets. 

The cases where the full pattern can be explicitly determined are clearly special. For example, one 
of the conditions used for exact characterization of patterns in this paper is that inside each patch 
the height variables are periodic and hence Ap (r) is constant. It is easy to check that this condition 
is not met for most of the sink geometries. For example, patterns of the type discussed in Section 
4 with any uj other than integer multiples of 7r/4 has aperiodic patches. Similarly for the patterns 
with two sources even slight deviation of the position of the second source in Fig llll from the x-axis 
introduces aperiodic patches. One of such patterns produced by adding 40000 grains each at (—180, 0) 
and (180, 20) is shown in FiglTTl The regions with stipes of red and yellow are the aperiodic patches. 
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Fig. 17 Pattern produced by adding — 40000 grains each at ( — 180, 0) and (180, 20) on F-lattice with initial 
chequerboard distribution of grains and relaxing. Color code red=0 and yellow=l. (Details can be seen in the 
online version using zoom in ) 



Also boundary of all these aperiodic patches have slopes other than 0, ±1 and oo and some boundaries 
of patches are not straight lines. In such cases, the present treatment is clearly not applicable. However 
the scaling analysis for the growth of spatial lengths in the pattern with N is still valid. 

The function D = d + ie satisfies discrete Cauchy-Riemann condition (equation ([35])). These func- 
tions are known as discrete holomorphic functions in the mathematics literature. They have been 
usually studied for a square grid of points on the plane [T8lfT9] . While more general discretizations of 
the plane have been discussed p01l21| . not much is known about the behavior of such functions for 
muti-sheeted Riemann surfaces. 

In our analysis we have also used a fact that the patterns have non-zero average overall excess 
density ( i.e. C2 in Eq. (8) is nonzero). The case when C2 is zero is quite different, and requires a 
substantially different treatment. We hope to discuss such patterns in a future publication P^ . 

Acknowledgements DD thanks Prof. A. Libchaber for suggesting this problem, and some useful discussions, 
and the Department of Science and Technology, India for financial support through a J.C. Bose fellowship. 
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